library(tidyverse)
library(plyr)
#install.packages("lfe")
library(lfe)
library(tidyr)
library(BBmisc)
library(dplyr)
library(broom)
library(corrplot)
library(RColorBrewer)
library(stargazer)
library(xtable)
library(mfx)
library(sjPlot)
library(texreg)
library(sjmisc)
library(interflex)


mytheme <- theme(
  panel.background = element_blank(),
  # panel.grid.major.y = element_line(color = "black"),
  strip.background = element_blank(),
  strip.text = element_text(face = "bold", size = 11),
  # text = element_text(family = "TT Times New Roman"),
  plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
  legend.key = element_blank(),
  axis.title.y = element_text(vjust = 3),
  axis.title.x = element_text(vjust = 0),
  panel.grid.major = element_line(color = "#cccccc", size = .2),
  axis.ticks = element_blank()
)

#=====================================================
#  Data Import
#=====================================================

cd<-"[[INSERT PATH HERE]]"

input<-str_c(cd)
graphs<-str_c(cd)
latex<-str_c(cd)


files<-list.files(path=input)

for (f in 1:length(files)){
  temp<-read.csv(paste(input,files[[f]],sep=""),encoding="UTF-8")
  temp<-temp%>%slice(3:nrow(temp))
  temp["db"]=f
  if (f==1){
    final<-temp
  }else{
    final<-bind_rows(final,temp)
  }
}

dta<-final%>%filter(Progress==100)

#attention check: only 10% failures, half of whom are people who did not leave the mail 
dta_good<-dta%>%filter(between(as.numeric(as.character(Q_Immigration_4)), 60, 70))

dta_good[,]<- lapply(dta_good, function(x) type.convert(as.character(x), as.is = TRUE))


#=====================================================
#  Create variables
#=====================================================
dta_good$Study<-as.numeric(revalue(dta_good$Study,c("Scuola superiore"=2,"Licenza media"=1,"Diploma universitario (triennale)"=3,"Diploma universitario (magistrale)"=4,"PhD / MBA o altri titoli post-universitari"=5)))
dta_good$Often_Offline<-as.numeric(revalue(dta_good$Often_Offline,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta_good$Often_Online<-as.numeric(revalue(dta_good$Often_Online,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta_good$Facebook<-as.numeric(revalue(dta_good$Facebook,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta_good$Internet<-as.numeric(revalue(dta_good$Internet,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta_good$Twitter<-as.numeric(revalue(dta_good$Twitter,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta_good$Like_GVT<-as.numeric(revalue(dta_good$Like_GVT,c("Completamente in disaccordo"=1,"In disaccordo"=2,"Abbastanza in disaccordo" =3,"Né d'accordo né in disaccordo"=4,"Abbastanza d'accordo"=5,"D'accordo" =6,"Completamente d'accordo"=7)))
dta_good$Pro_Gov<-ifelse(dta_good$Like_GVT==4,0,ifelse(dta_good$Like_GVT>4,1,-1))
dta_good$Age<-as.numeric(dta_good$Age)

dta_good<-dta_good%>%
  filter(Age>17)

dta_good$Like_GVT_<-ifelse(dta_good$Like_GVT==1,"Strongly Disagree",
                           ifelse(dta_good$Like_GVT==2,"Disagree",
                                  ifelse(dta_good$Like_GVT==3,"Tend to Disagree",
                                         ifelse(dta_good$Like_GVT==4,"Neutral",
                                                ifelse(dta_good$Like_GVT==5,"Tend to Agree",
                                                       ifelse(dta_good$Like_GVT==6,"Agree",
                                                              ifelse(dta_good$Like_GVT==7,"Strongly Agree",NA)))))))


dta_good$Pro_Gov_<-ifelse(dta_good$Pro_Gov==-1,"Anti Government",
                          ifelse(dta_good$Pro_Gov==0,"Neutral","Pro Government"))
dta_good$Sex<-0
dta_good$Sex[dta_good$Gender=="Uomo"]<-1

#number for all tech savvy and compute overall savvyness
dta_good$Savvy<-0
for (i in grep("Tech_Savv",names(dta_good))){
  assign("var",names(dta_good)[[i]])
  dta_good[[var]]<-as.numeric(revalue(dta_good[[var]],c("1: Non conosco"=1,"2"=2,"3"=3,"4"=4,"5: Conosco bene"=5)))
  dta_good$Savvy<-dta_good$Savvy+as.numeric(dta_good[[var]])
}
dta_good$Savvy<-dta_good$Savvy/(13*5)

#=====================================================
#  Create variables for whole dataset 
#=====================================================

library(plyr)
dta$Study<-as.numeric(revalue(dta$Study,c("Scuola superiore"=2,"Licenza media"=1,"Diploma universitario (triennale)"=3,"Diploma universitario (magistrale)"=4,"PhD / MBA o altri titoli post-universitari"=5)))
dta$Often_Offline<-as.numeric(revalue(dta$Often_Offline,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta$Often_Online<-as.numeric(revalue(dta$Often_Online,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta$Facebook<-as.numeric(revalue(dta$Facebook,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta$Internet<-as.numeric(revalue(dta$Internet,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta$Twitter<-as.numeric(revalue(dta$Twitter,c("Meno spesso"=1,"Qualche volta al mese"=2,"1 o 2 volte alla settimana"=3,"3 - 6 volte alla settimana"=4,"Una volta al giorno"=5,"Più volte al giorno"=6,"Praticamente sempre"=7)))
dta$Like_GVT<-as.numeric(revalue(dta$Like_GVT,c("Completamente in disaccordo"=1,"In disaccordo"=2,"Abbastanza in disaccordo" =3,"Né d'accordo né in disaccordo"=4,"Abbastanza d'accordo"=5,"D'accordo" =6,"Completamente d'accordo"=7)))
dta$Pro_Gov<-ifelse(dta$Like_GVT==4,0,ifelse(dta$Like_GVT>4,1,-1))
dta$Age<-as.numeric(dta$Age)

dta$Sex<-0
dta$Sex[dta$Gender=="Uomo"]<-1

#number for all tech savvy and compute overall savvyness
dta$Savvy<-0
for (i in grep("Tech_Savv",names(dta))){
  assign("var",names(dta)[[i]])
  dta[[var]]<-as.numeric(revalue(dta[[var]],c("1: Non conosco"=1,"2"=2,"3"=3,"4"=4,"5: Conosco bene"=5)))
  dta$Savvy<-dta$Savvy+as.numeric(dta[[var]])
}
dta$Savvy<-dta$Savvy/(13*5)

dta$attention<-ifelse(between(as.numeric(dta$Q_Immigration_4),60,70),1,0)

detach(plyr)

dta_good$age_reg=dta_good$age_

#=====================================================
#  Create Variables CB score
#=====================================================

##work on CB choice pre-treatment
#as.numeric choice rank
#drop one observation with missing choice for CB (weird)
dta_good<-dta_good%>%filter(Choice_2_1!='')

dta_good$Choice_1_2=as.numeric(as.character(dta_good$Choice_1_2))
dta_good$Choice_1_3=as.numeric(as.character(dta_good$Choice_1_3))
dta_good$Choice_1_4=as.numeric(as.character(dta_good$Choice_1_4))
dta_good$Choice_1_5=as.numeric(as.character(dta_good$Choice_1_5))
dta_good$Choice_1_8=as.numeric(as.character(dta_good$Choice_1_8))

dta_good$Choice_2_1_p=as.numeric(as.character(dta_good$Choice_2_1))
dta_good$Choice_2_2_p=as.numeric(as.character(dta_good$Choice_2_2))
dta_good$Choice_2_3_p=as.numeric(as.character(dta_good$Choice_2_3))
dta_good$Choice_2_4_p=as.numeric(as.character(dta_good$Choice_2_4))
dta_good$Choice_2_5_p=as.numeric(as.character(dta_good$Choice_2_5))

#questions are always in this order CB-PRO CB-AG NCB-PRO NCB-AG NEU
dta_good["CB_reader_N"]=0

dta_good$CB_reader_N[dta_good$Choice_1_2==1]=dta_good$CB_reader_N[dta_good$Choice_1_2==1]+1
dta_good$CB_reader_N[dta_good$Choice_1_3==1]=dta_good$CB_reader_N[dta_good$Choice_1_3==1]+1
dta_good$CB_reader_N[dta_good$Choice_2_1==1]=dta_good$CB_reader_N[dta_good$Choice_2_1==1]+1
dta_good$CB_reader_N[dta_good$Choice_2_2==1]=dta_good$CB_reader_N[dta_good$Choice_2_2==1]+1

dta_good <- dta_good%>%
  mutate(
    CF_reader_once = ifelse(CB_reader_N>0,1,0)
  )

#create score --> highest score = ranked higher news of type k
#category k gets 5 points if the first news was in that group, 4 for 2nd....
dta_good["CB_reader_Score"]=20-dta_good$Choice_1_2-dta_good$Choice_1_3-dta_good$Choice_2_1-dta_good$Choice_2_2

#=====================================================
#  Analysis incentive treatment (chosen)
#=====================================================
#=====================================================
#  Create Variables
#=====================================================

#create treatment for incentive in choosing article
dta_good["treatment_chosen"]=-1
dta_good$treatment_chosen[dta_good$TR_LOWINC=="Ho capito il funzionamento della prossima sezione."]=1
dta_good$treatment_chosen[dta_good$TR_HIGHINC=="Ok, ho capito il funzionamento della prossima sezione."]=2
dta_good$treatment_chosen[dta_good$TR_NOINC=="Ok, vai alla prossima sezione."]=0

#variable for chosen article
#CB are more often chosen than NCB
dta_good["chosen"]=-1
dta_good$chosen[dta_good$Choice_Chosen_1==1]=1
dta_good$chosen[dta_good$Choice_Chosen_2==1]=2
dta_good$chosen[dta_good$Choice_Chosen_3==1]=3
dta_good$chosen[dta_good$Choice_Chosen_4==1]=4
dta_good$chosen[dta_good$Choice_Chosen_5==1]=5

#absolutely no variation across treatments
table(dta_good$chosen,dta_good$treatment_chosen)

dta_good["choseCB"]=0
dta_good$choseCB[dta_good$chosen==1]=1
dta_good$choseCB[dta_good$chosen==2]=1
dta_good["choseNCB"]=0
dta_good$choseNCB[dta_good$chosen==3]=1
dta_good$choseNCB[dta_good$chosen==4]=1
dta_good["choseNEU"]=0
dta_good$choseNEU[dta_good$chosen==5]=1
dta_good["choseCBPRO"]=0
dta_good$choseCBPRO[dta_good$chosen==1]=1
dta_good["choseAG"]=0
dta_good$choseAG[dta_good$chosen==2]=1
dta_good$choseCB[dta_good$chosen==4]=1
dta_good["chosePRO"]=0
dta_good$chosePRO[dta_good$chosen==1]=1
dta_good$chosePRO[dta_good$chosen==3]=1
dta_good["choseNONE"]=0
dta_good$choseNONE[dta_good$chosen==-1]=1

dta_good$choseSAME<-ifelse((dta_good$Pro_Gov==1&dta_good$chosePRO==1)|(dta_good$Pro_Gov==-1&dta_good$choseAG==1),1,
                           ifelse((dta_good$Pro_Gov==1&dta_good$choseAG==1)|(dta_good$Pro_Gov==-1&dta_good$chosePRO==1),0,NA))

dta_good["age_"]=NA
dta_good$age_[as.numeric(dta_good$Age)<20]="young"
dta_good$age_[between(as.numeric(dta_good$Age),20,40)]="young_adult"
dta_good$age_[between(as.numeric(dta_good$Age),41,60)]="adult"
dta_good$age_[as.numeric(dta_good$Age)>60]="old"

dta_good$educ_=ifelse(dta_good$Study==1,"Low",
                      ifelse(dta_good$Study==2,"Medium", "High"))

dta_good$no_incentive<-ifelse(dta_good$treatment_chosen==0,1,0)

dta_good$TIME_ARTICLE_CHOSEN_Page.Submit=as.numeric(dta_good$TIME_ARTICLE_CHOSEN_Page.Submit)

dta_good$treatment_chosen_label=ifelse(dta_good$treatment_chosen==0,"No Inc",
                                       ifelse(dta_good$treatment_chosen==1,"Low Inc",
                                              ifelse(dta_good$treatment_chosen==2,"High Inc",NA)))

dta_good$treatment_chosen_label_2 = factor(dta_good$treatment_chosen_label,levels=c("No Inc","Low Inc","High Inc")) #plot with correct labels
dta_good$educ_level = factor(dta_good$educ_,levels=c("Low","Medium","High"))

perc.rank.savvy <- ecdf(dta_good$Savvy)
dta_good$Savvy_perc=perc.rank.savvy(dta_good$Savvy)*100
dta_good$Savvy_<-ifelse(dta_good$Savvy_perc<25,"Bottom25",
                        ifelse(dta_good$Savvy_perc<75,"Middle50","Top25"))

#=========================================================================
#  Export selected figures for text
#=========================================================================

controls=c("CB_reader_Score" ,"Educ","Age","Like_GVT","Sex","Internet","Facebook","Twitter","Often_Online","Often_Offline","Savvy","log(TIME_ARTICLE_CHOSEN_Page.Submit)")
controls_light=c("CB_reader_Score" ,"Educ","Age","Sex","Internet","Facebook","Twitter","Often_Online","Often_Offline","Savvy","log(TIME_ARTICLE_CHOSEN_Page.Submit)")

#======
# chosen: choice
# Descriptive_Chosen
#======
stats_to_export<-dta_good%>%
  dplyr::group_by(chosen)%>%
  dplyr::summarise(
    Men=mean(Sex,na.rm=TRUE),
    Age=mean(Age,na.rm=TRUE),
    Education=mean(Study,na.rm=TRUE),
    Like_GVT=mean(Like_GVT,na.rm=TRUE),
    Internet=mean(Internet,na.rm=TRUE),
    Facebook=mean(Facebook,na.rm=TRUE),
    Twitter=mean(Twitter,na.rm=TRUE),
    News_Offline=mean(Often_Offline,na.rm=TRUE),
    News_Online=mean(Often_Online,na.rm=TRUE),
    Tech_Savvy= mean(Savvy,na.rm=TRUE),
    n_obs=n(),
    Chosen_Article=first(chosen)
  )%>%
  dplyr::select(Chosen_Article,Men,Age,Education,Like_GVT,Internet,Facebook,Twitter,News_Offline,News_Online,Tech_Savvy,n_obs)%>%
  mutate(
    Chosen_Article=ifelse(Chosen_Article==-1,"None",
                          ifelse(Chosen_Article==1,"CB-Pro",
                                 ifelse(Chosen_Article==2,"CB-Against",
                                        ifelse(Chosen_Article==3,"Non CB-Pro",
                                               ifelse(Chosen_Article==4,"Non CB-Against","Neutral")))))
  )

stats_to_export_allpop<-dta_good%>%
  dplyr::group_by(Progress)%>%
  dplyr::summarise(
    Men=mean(Sex,na.rm=TRUE),
    Age=mean(Age,na.rm=TRUE),
    Education=mean(Study,na.rm=TRUE),
    Like_GVT=mean(Like_GVT,na.rm=TRUE),
    Internet=mean(Internet,na.rm=TRUE),
    Facebook=mean(Facebook,na.rm=TRUE),
    Twitter=mean(Twitter,na.rm=TRUE),
    News_Offline=mean(Often_Offline,na.rm=TRUE),
    News_Online=mean(Often_Online,na.rm=TRUE),
    Tech_Savvy= mean(Savvy,na.rm=TRUE),
    n_obs=n(),
    Chosen_Article="All"
  )%>%
  dplyr::select(Chosen_Article,Men,Age,Education,Like_GVT,Internet,Facebook,Twitter,News_Offline,News_Online,Tech_Savvy,n_obs)

stats_to_export=rbind(stats_to_export,stats_to_export_allpop)

#transpose table to export
n <- stats_to_export$Chosen_Article
# transpose all but the first column (name)
stats_to_export <- as.data.frame(t(stats_to_export[,-1]))
colnames(stats_to_export) <- n
#stats_to_export$myfactor <- factor(row.names(stats_to_export))
#====
# Table 1
#====


print(xtable(stats_to_export, type = "latex"), file = str_c(latex,"/Descriptive_Chosen.tex"))

#====
# CB who chooses it?
# CBScore_who
#====

myols.who <- lm(CB_reader_Score ~Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                data = dta_good)
summary(myols.who)

probit.who <- probitmfx(CF_reader_once ~ Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                        data = dta_good)

writeLines(capture.output(stargazer(myols.who,probit.who, title="Taste for CB", align=TRUE,
                                    dep.var.labels=c("OLS","Probit"),
                                    covariate.labels=c("Educ","Age","Sex","Internet","Facebook","Twitter","LIke GVT","Often Online","Often Offline","Savvy"))), str_c(latex,"/CBScore_who.tex"))


#====
# Table 2
#====

writeLines(capture.output(texreg(
  list(myols.who,probit.who),
  custom.coef.names = c("Constant","Education","Age","Sex","Internet","Facebook","Twitter","Like GVT","Often Online","Often Offline","Dig Literacy"),
  caption = "Taste for CB",
  #custom.header = c(""),
  custom.model.names = c("OLS","Probit"),
  label = "CBScore_who",
  stars = c(0.1, 0.05, 0.01),
  digits =3,
  dcolumn = TRUE,
  booktabs = TRUE,
  use.packages = FALSE,
  include.adjrs = FALSE,
  include.bic = FALSE,
  include.aic = FALSE,
  str_c(latex,"/CBScore_who.tex"))))


#====
# Chosen: time spent
# Effort_Chosen_Read
#====

dta_good$wordcount_chosen=ifelse(dta_good$chosen==1,193,
                                 ifelse(dta_good$chosen==2,187,
                                        ifelse(dta_good$chosen==3,160,
                                               ifelse(dta_good$chosen==4,161,
                                                      ifelse(dta_good$chosen==5,157,NA)))))

myols.notreat <- lm(log(TIME_CHOICE_CHOSEN_Page.Submit+1) ~ CB_reader_Score +Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                    data = dta_good)

myols.nocontrols <- lm(log(TIME_CHOICE_CHOSEN_Page.Submit+1) ~ as.factor(treatment_chosen), 
                       data = dta_good)

myols.treat <- lm(log(TIME_CHOICE_CHOSEN_Page.Submit+1) ~ as.factor(treatment_chosen) +CB_reader_Score +Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                  data = dta_good)

myols.readnotreat <- lm(log(TIME_ARTICLE_CHOSEN_Page.Submit/wordcount_chosen+1) ~ CB_reader_Score +Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                        data = dta_good)

myols.readnocontrols <- lm(log(TIME_ARTICLE_CHOSEN_Page.Submit/wordcount_chosen+1) ~ as.factor(treatment_chosen), 
                           data = dta_good)

myols.readtreat <- lm(log(TIME_ARTICLE_CHOSEN_Page.Submit/wordcount_chosen+1) ~ as.factor(treatment_chosen) +CB_reader_Score +Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                      data = dta_good)

#====
# Table 3
#====

writeLines(capture.output(stargazer(myols.notreat,myols.treat,myols.nocontrols,myols.readnotreat,myols.readtreat,myols.readnocontrols, title="Effort under different treatments", align=TRUE,
                                    dep.var.labels=c("Time spent choosing the article","Time spent reading the chosen article"),
                                    covariate.labels=c("Low Inc","High Inc","CB Score","Educ","Age","Sex","Internet","Facebook","Twitter","LIke GVT","Often Online","Often Offline","Savvy"))), str_c(latex,"/Effort_Chosen_Read.tex"))


#############

summary(dta_good$TIME_ARTICLE_CHOSEN_Page.Submit)

dta_good$treatment_chosen_label=ifelse(dta_good$treatment_chosen==0,"No Inc",
                                       ifelse(dta_good$treatment_chosen==1,"Low Inc",
                                              ifelse(dta_good$treatment_chosen==2,"High Inc",NA)))


#====
# Chosen: choseCB treatment
# Likelihood_CB_combined
#====

baseline <- probitmfx(choseCB ~ as.factor(treatment_chosen) +CB_reader_Score +Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                      data = dta_good)
myprobit.Treat_age <- probitmfx(choseCB ~ as.factor(treatment_chosen)*Age +CB_reader_Score +Study+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                                data = dta_good)
myprobit.Treat_edubins <- probitmfx(choseCB ~ as.factor(treatment_chosen)*educ_ +CB_reader_Score +Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                                    data = dta_good)
myprobit.Treat_savbins <- probitmfx(choseCB ~ as.factor(treatment_chosen)*Savvy_ +Age +CB_reader_Score+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, 
                                    data = dta_good)
myprobit.all_together <- probitmfx(choseCB ~ as.factor(treatment_chosen)*Age +as.factor(treatment_chosen)*educ_ +as.factor(treatment_chosen)*Savvy_ +Age +CB_reader_Score+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, 
                                   data = dta_good)
#====
# Table 4
#====

writeLines(capture.output(texreg(
  list(baseline,myprobit.Treat_age,myprobit.Treat_edubins,myprobit.Treat_savbins,myprobit.all_together),
  #custom.coef.names = c("CB Score","Educ","Age","Sex","Internet","Facebook","Twitter","LIke GVT","Often Online","Often Offline","Savvy","Low Inc","High Inc"),
  caption = "Likelihood of choosing CB",
  ##custom.header = c("CB"),
  custom.model.names = c("(1)","(2)","(3)","(4)","(5)"),
  label = "Likelihood_CB_combined",
  stars = c(0.1, 0.05, 0.01),
  digits =3,
  dcolumn = TRUE,
  booktabs = TRUE,
  use.packages = FALSE,
  include.adjrs = FALSE,
  include.bic = FALSE,
  include.aic = FALSE,
  custom.note = ("\\parbox{.8\\linewidth}{\\vspace{2pt}%stars. \\\\
                 Probit model on the Likelihood of choosing a CB article\\Marginal effects and standard error in brackets}"))), str_c(latex,"/Likelihood_CB_combinedTables.tex"))


#====
# Chosen: choseCB treatment
# Likelihood_CB_combined_percentile
#====

perc.rank.savvy <- ecdf(dta_good$Savvy)
dta_good$savvy_perc=perc.rank.savvy(dta_good$Savvy)*100

perc.rank.age <- ecdf(dta_good$Age)
dta_good$age_perc=perc.rank.age(dta_good$Age)*100

perc.rank.educ <- ecdf(dta_good$Study)
dta_good$educ_perc=perc.rank.educ(dta_good$Study)*100


baseline <- probitmfx(choseCB ~ as.factor(treatment_chosen) +CB_reader_Score +educ_perc+age_perc+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+savvy_perc, 
                      data = dta_good)
myprobit.Treat_age <- probitmfx(choseCB ~ as.factor(treatment_chosen)*age_perc +CB_reader_Score +educ_perc+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+savvy_perc, 
                                data = dta_good)
myprobit.Treat_edubins <- probitmfx(choseCB ~ as.factor(treatment_chosen)*educ_perc +CB_reader_Score +age_perc+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+savvy_perc, 
                                    data = dta_good)
myprobit.Treat_savbins <- probitmfx(choseCB ~ as.factor(treatment_chosen)*savvy_perc +age_perc +CB_reader_Score+educ_perc+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, 
                                    data = dta_good)
myprobit.all_together <- probitmfx(choseCB ~ as.factor(treatment_chosen)*age_perc +as.factor(treatment_chosen)*educ_perc +as.factor(treatment_chosen)*savvy_perc  +CB_reader_Score+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, 
                                   data = dta_good)

#====
# Figure 1
#====

writeLines(capture.output(texreg(
  list(baseline,myprobit.Treat_age,myprobit.Treat_edubins,myprobit.Treat_savbins,myprobit.all_together),
  #custom.coef.names = c("CB Score","Educ","Age","Sex","Internet","Facebook","Twitter","LIke GVT","Often Online","Often Offline","Savvy","Low Inc","High Inc"),
  caption = "Likelihood of choosing CB",
  #custom.header = c("CB"),
  custom.model.names = c("(1)","(2)","(3)","(4)","(5)"),
  label = "Likelihood_CB_combined_percentile",
  stars = c(0.1, 0.05, 0.01),
  digits =3,
  dcolumn = TRUE,
  booktabs = TRUE,
  use.packages = FALSE,
  include.adjrs = FALSE,
  include.bic = FALSE,
  include.aic = FALSE,
  custom.note = ("\\parbox{.8\\linewidth}{\\vspace{2pt}%stars. \\\\
                 Probit model on the Likelihood of choosing a CB article\\\\Marginal effects and standard error in brackets}"))), str_c(latex,"/Likelihood_CB_combinedTables_percentile.tex"))


#====
# Chosen: choseSAME treatment
# Likelihood_SAME_combined
#====


baseline <- probitmfx(choseSAME ~ as.factor(treatment_chosen) +CB_reader_Score +Study+Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                      data = dta_good)
myprobit.Treat_age <- probitmfx(choseSAME ~ as.factor(treatment_chosen)*Age +CB_reader_Score +Study+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                                data = dta_good)
myprobit.Treat_edubins <- probitmfx(choseSAME ~ as.factor(treatment_chosen)*educ_ +CB_reader_Score +Age+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, 
                                    data = dta_good)
myprobit.Treat_savbins <- probitmfx(choseSAME ~ as.factor(treatment_chosen)*Savvy_ +Age +CB_reader_Score+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, 
                                    data = dta_good)
myprobit.all_together <- probitmfx(choseSAME ~ as.factor(treatment_chosen)*Age +as.factor(treatment_chosen)*educ_ +as.factor(treatment_chosen)*Savvy_ +Age +CB_reader_Score+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, 
                                   data = dta_good)
#====
#Table 6
#====
writeLines(capture.output(texreg(
  list(baseline,myprobit.Treat_age,myprobit.Treat_edubins,myprobit.Treat_savbins,myprobit.all_together),
  #custom.coef.names = c("CB Score","Educ","Age","Sex","Internet","Facebook","Twitter","LIke GVT","Often Online","Often Offline","Savvy","Low Inc","High Inc"),
  caption = "Likelihood of choosing CB",
  #custom.header = c("CB"),
  custom.model.names = c("(1)","(2)","(3)","(4)","(5)"),
  label = "Likelihood_CB_combined",
  stars = c(0.1, 0.05, 0.01),
  digits =3,
  dcolumn = TRUE,
  booktabs = TRUE,
  use.packages = FALSE,
  include.adjrs = FALSE,
  include.bic = FALSE,
  include.aic = FALSE,
  custom.note = ("\\parbox{.8\\linewidth}{\\vspace{2pt}%stars. \\\\
                 Probit model on the Likelihood of choosing a CB article\\Marginal effects and standard error in brackets}"))), str_c(latex,"/Likelihood_SAME_combinedTables.tex"))


#=========================================================
# Figures
#=========================================================

myprobit.Treat_age <- glm(choseCB ~ as.factor(treatment_chosen_label)*poly(Age,3) +CB_reader_Score +Study+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline+Savvy, family = binomial(link = "probit"), 
                          data = dta_good[!is.na(dta_good$Age)&between(dta_good$Age,18,70)&dta_good$treatment_chosen_label!="Low Inc",])

#====
# Figure 1 
#====
plot_model(myprobit.Treat_age, type = "pred", terms = c("Age","treatment_chosen_label")) +theme_light()+
  ggtitle("Predicted probabilities of choosing CB") + labs(color='Treatment') +
  ylab("Probability of choosing CB")
ggsave(paste(graphs,"/Interaction_TreatAge_Poly3_noLow.png",sep=''),plot=last_plot())


myprobit.Treat_educ <- glm(choseCB ~ as.factor(treatment_chosen_label_2)*educ_level+Savvy+Age +CB_reader_Score +Like_GVT+Sex+Internet+Facebook+Twitter+Often_Online+Often_Offline, family = binomial(link = "probit"), 
                           data = dta_good[!is.na(dta_good$treatment_chosen_label_2)&dta_good$treatment_chosen_label_2!="Low Inc",])
plot_model(myprobit.Treat_educ, type = "pred", terms = c("educ_level","treatment_chosen_label_2")) +theme_light()+
  ggtitle("Predicted probabilities of choosing CB") + labs(color='Treatment') +
  ylab("Probability of choosing CB") + xlab("Level of Education")
ggsave(paste(graphs,"/Interaction_TreatEduc.png",sep=''),plot=last_plot())

myprobit.Treat_savvy <- glm(choseCB ~ as.factor(treatment_chosen_label)*poly(Savvy_perc,3)+Age +CB_reader_Score +Study+Study+Sex+Internet+Facebook+Twitter+Like_GVT+Often_Online+Often_Offline, family = binomial(link = "probit"), 
                            data = dta_good[!is.na(dta_good$Savvy_perc)&dta_good$treatment_chosen_label!="Low Inc",])
plot_model(myprobit.Treat_savvy, type = "pred", terms = c("Savvy_perc","treatment_chosen_label")) +theme_light()+
  ggtitle("Predicted probabilities of choosing CB") + labs(color='Treatment') +
  ylab("Probability of choosing CB") + xlab("Percentile of Digital Literacy")
ggsave(paste(graphs,"/Interaction_TreatSavvy_Poly3_noLow.png",sep=''),plot=last_plot())

#=================
# Interflex
#=================

dta_good <- dta_good%>%
  mutate(
    treat_bin = case_when(treatment_chosen_label == "High Inc" ~ 1,
                          treatment_chosen_label == "No Inc" ~ 0,
                          TRUE ~ NA_real_)
  )


#====
# Figure 3
#====
inter.kernel(Y="choseCB",D="treat_bin",X="Age",Ylabel = "Clickbait",Z=c("CB_reader_Score","Study","Sex","Internet","Facebook","Twitter","Like_GVT","Often_Online","Often_Offline","Savvy"),na.rm = TRUE, Dlabel = "Monetary Incentives", Xlabel="Age", theme.bw = TRUE,main = "",data=dta_good[!is.na(dta_good$Age)&between(dta_good$Age,18,70)&dta_good$treatment_chosen_label!="Low Inc",])
ggsave(paste(graphs,"/Interkernel_AGE_ALLcontrols.png",sep=''),plot=last_plot())

inter.kernel(Y="choseCB",D="treat_bin",X="Study",Ylabel = "Clickbait",Z=c("CB_reader_Score","Age","Sex","Internet","Facebook","Twitter","Like_GVT","Often_Online","Often_Offline","Savvy"),na.rm = TRUE, Dlabel = "Monetary Incentives", Xlabel="Education", theme.bw = TRUE,main = "",data=dta_good[!is.na(dta_good$Age)&between(dta_good$Age,18,70)&dta_good$treatment_chosen_label!="Low Inc",])
ggsave(paste(graphs,"/Interkernel_EDUC_ALLcontrols.png",sep=''),plot=last_plot())

inter.kernel(Y="choseCB",D="treat_bin",X="Savvy",Ylabel = "Clickbait",Z=c("CB_reader_Score","Age","Sex","Internet","Facebook","Twitter","Like_GVT","Often_Online","Often_Offline","Study"),na.rm = TRUE, Dlabel = "Monetary Incentives", Xlabel="Digital Literacy", theme.bw = TRUE,main = "",data=dta_good[!is.na(dta_good$Age)&between(dta_good$Age,18,70)&dta_good$treatment_chosen_label!="Low Inc",])
ggsave(paste(graphs,"/Interkernel_DigLitPercentile_ALLcontrols.png",sep=''),plot=last_plot())